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We report an extension of the source imaging method for analyzing three-dimensional sources 
from three-dimensional correlations. Our technique consists of expanding the correlation data and 
the underlying source function in spherical harmonics and inverting the resulting system of one- 
dimensional integral equations. With this strategy, we can image the source function quickly, even 
with the finely binned data sets common in three-dimensional analyses. 

PACS numbers: PACS numbers: 25.75.-q, 25.75.Gz 



INTRODUCTION 

The Bose-Einstein induced correlation of particles pro- 
duced in high-energy pp reactions was first observed by 
Goldhaber, Goldhaber, Lee and Pais Q in 1960. The 
connection of this effect to Hanbury-Brown/Twiss in- 
tensity interferometry followed several years later 0, Q • 
Since then, two particle interferometry has become a 
powerful tool to study heavy-ion reactions: it is used 
at all available energies from MSU-NSCL 3 to RH1C 
[j| and the Tevatron [||. The variety of particles used 
in these studies include protons, neutrons, mesons, in- 
termediate mass fragments, electrons, photons, and even 
exotic hadrons such as A and S's. In all cases, our goal is 
to understand the space-time development of nuclear re- 
actions. One particularly fruitful class of this observable 
is multidimensional pion correlations. Two pion correla- 
tions have been used to demonstrate the hydrodynamical 
scaling of the source radii 0,0, to determine the average 
freeze-out phase-space density 0, 0, and to charac- 
terize the final state of heavy-ion collisions 0, Q. We 
would like to extend this success to other particles, such 
as protons or kaons, where simple Coulomb corrections 
may be (or are) insufficient to study the data. Indeed, 
final-state interactions (FSI) are essential ingredients 0| 
for multidimensional measurements of non-identical cor- 
relations. For this, we extend the one dimensional source 
imaging technique of 0, [Til fl5| to these multidimen- 
sional correlations. 

Source imaging allows us to cleanly separate effects due 
to FSI and symmetrization from effects due to the source 
function itself. Because no assumptions are made about 
the functional form of the source function, the imaging 
is model independent. This means that it is possible to 
image non-Gaussian sources jl6j . Furthermore, imaging 
allows us to directly compare full three-dimensional pro- 
ton and pion sources, free from the distortions due to 
FSI. Imaging is a useful alternate approach to the stan- 
dard analysis which involves first applying a Coulomb 
correction then fitting the corrected data to a Gaussian. 



We can extract the source without an ad-hoc Coulomb 
correction, and later if we wish, fit the extracted source 
to a Gaussian. While this technique may be applied to 
other pair types, e.g. proton-proton pairs |l7j . we will 
concentrate on like meson pairs. 

Extracting the source function, Sp(r'), for a pair of 
identical particles begins by noting that the Koonin- 
Pratt equation relates it to the experimentally measured 
two-particle correlation, Cp(q') (or T^p(q'), its deviation 
from 1) [30: 



ftp(q') = Cp(q') -l=f dr'A'(q',r') Sp(r') 



(1) 



Thus, "imaging the source" means inverting this equa- 
tion. In JJJ, P = pi + P2 is the total momentum of the 
pair and the kernel of the integral equation is 



A-(q',r') = |<7 ) (r')| 2 -l 



(2) 



The wavefunction, describes the propagation of the 

pair from a relative separation of r' in the pair center of 
mass (CM) to the detector with relative momentum q' = 
^(p'i — Pa)- Here primes denote quantities in the pair 
CM frame. The source function itself is the probability 
of emitting the pair a distance of r' apart in the pair CM 
frame. 

The problem of inverting Eq. is generally ill-posed, 
meaning that small fluctuations in the data, even if well 
within statistical or systematic errors, can lead to large 
changes in the imaged source function. In angle-averaged 
correlations, this stability problem is solved by the use 
of constraints, a basis spline representation of the source, 
and "optimized discretization" 0, Q, . By expand- 
ing the source and the correlation in spherical harmonics, 
we may convert the full three-dimensional problem into a 
series of one-dimensional problems where these tools are 
equally applicable. Throughout this process, we must 
take care to respect the frame dependence of the source 
function. In an angle-averaged like-pair correlation, one 
may work in terms of qw = \J q 2 — Qq , so the analysis 
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frame is not an issue. In full three-dimensional analysis, 
we must properly carry out the transformation from the 
frame where the data is compiled to the frame where we 
image. 

We now outline this paper. First, we will set up the 
problem, taking into account the special issues in a full 
three-dimensional analysis. This includes a discussion of 
how the source function and underlying emission func- 
tions are related, the frame dependence of the imaging, 
the coordinate system used in the imaging and a non- 
trivial test case. Next, we describe the three-dimensional 
imaging procedure. This process amounts to breaking 
the problem into a series of one-dimensional problems. 
We describe this procedure and its implications, espe- 
cially to the source and its representation. We illustrate 
the entire imaging process using the non-trivial test case. 



PRELIMINARIES 

In this section, we detail the connection between the 
correlation and source functions and the underlying emis- 
sion functions. In this process, we illuminate some of the 
technical issues arising from the frame dependence of the 
source function. We also define the coordinates that we 
use in our analyses. Following this, we detail the non- 
trivial model correlation we use as a test case to illustrate 
our analysis procedure. 



Relating the Source Function and Emission 
Functions 



Before we can image, we write the Koonin-Pratt equa- 
tion in a covariant form. This is important because the 
imaging is simplest in the pair CM frame, yet the data is 
often tabulated in an entirely different frame. We begin 
with the measured correlation which written as a con- 
volution of the final state wavefunction with a relative 
distance distribution: 



C P (q) =E 1 E 2 



dpi dp 2 



E 



dN x 
dpi 



E-2 



dN 2 
dp 2 



d 4 r 



(3) 

Here the relative distance distribution, Dp(r), is the 
probability of emitting a pair with a space-time separa- 
tion r = (ro, r). Each particle has approximately momen- 
tum P/2 (in the smoothness approximation, the slight 
difference between using pi, P2 and P/2 is ignored ja]). 
The relative distance distribution can be written as the 
convolution of emission functions for particles 1 and 2: 

£> P (r)= J d A RD 1 (R+r/2,-p/2)D 2 (R-r/2,-p/2). (4) 



The emission functions are the normalized particle emis- 
sion rates: 



D(r,p) 



Ed 7 N I Ed 3 N 



d 4 r d 3 p I d 3 p 



(5) 



This rate could in principal be computed in models such 
as RQMD or the Blast- Wave 

To make contact with Eq. QJ, we boost to the pair 
CM frame. Since Eq. ® is written covariantly, we may 
re-express it in any frame simply by replacing coordi- 
nates in the original frame with coordinates in the pair 
CM frame. Further, since the the wavefunction is inde- 
pendent of the relative time, r' Q , in the pair CM frame 
[42| , the time integral may be brought over to act on the 
relative distance distribution: 



C P (q) = J dV |<Mr')| 2 | dr' V P (r) 



(0) 



This allows us to define the source function in the pair 
CM frame: 



Sp(r') = J dr' V P (r). 



(7) 



With this source, we arrive at the Koonin-Pratt equation 
in Eq. (Q|. 

The Coordinate System 

We use the Bertsch-Pratt coordinates 0, in both 
r and q spaces and we define these coordinates in the 
rest frame of the interacting nuclear system. Once the 
directions of the unit vectors are defined, we can use the 
coordinate system in any frame boosted relative to the 
system frame. In the Bertsch-Pratt coordinates, the lon- 
gitudinal direction = z is along the beam direction 
and the sideward direction eg is perpendicular to the 
longitudinal direction and the direction of the total pair 
momentum, i.e. §5 = P/|P| x ej, = x. The final direc- 
tion, the outward direction, lies in the plane defined by 
the beam direction and P. The unit vector in the out- 
ward direction &o may be obtained from the other two 
unit vectors: = ej, x % = y. Together, the three 
unit vectors define the right-handed coordinate system 
shown in Fig. ^ ln this figure, we also define the po- 
lar angles (9, <p) which we use when working in spherical 
coordinates. 

Our use of these coordinates does not preclude per- 
forming an azimuthal HBT-type analysis. In an az- 
imuthal analysis, the part of the total momentum vector 
P perpendicular to the longitudinal axis is further decom- 
posed into components parallel to and perpendicular to 
the reaction plane. Since the Bertsch-Pratt coordinates 
are defined relative to the P direction, a correlation (or 
source) in these coordinates can be measured (or imaged) 
with respect to the reaction plane. 
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FIG. 1: The Bertsch-Pratt coordinate system in q-space. 

Simple Test Problem 

To test the imaging codes, we need to generate a test 
source function and correlation. To do this, we first 
choose the single particle source D(r, t, p) to use in Eqs. 
JTJ and 0. Then we perform the integrals in Eq. 
using the Correlation AfterBurner (CRAB) code |23| to 
generate the correlations. We use the symmetrized rela- 
tive wavefunction, including Coulomb repulsion. 

In the lab frame, our single particle source is composed 
of a prompt source and a long-lived source: 

D(r, t, p) <x (/*(*) + (1 - f)9(t) exp (-t/r)) 

/ p 2 \ (8) 

Here, the Gaussian given by 

^aus S (r)ocexp(-^-^-^). (9) 

Our test single particle source is radially symmetric with 
Rx = Ry — Rz — 4 fm, with a life-time of r = 20 fm/c 
(comparable to the uj lifetime of 23 fm/c and the sys- 
tem lifetime in other models H^)- The fraction of pions 
emitted from the prompt versus long-lived source is con- 
trolled by / and is set to / = 0.65, corresponding roughly 
to the number of pions created from all resonance decays 
in reaction simulations [20I |25| . We set the source tem- 
perature to 50 MeV. 

To produce the source function in Fig. [5J we first gen- 
erated a set of phase-space points sampled from the sin- 
gle particle source in Eq. (JHJ). We then perform the 
integrals in Eqs. J3J and Q using this set of phase-space 
points. Since we did not have pairs with exactly the same 
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FIG. 2: The model source function constructed from the 
model emission functions in Eq. JHJ. 

momentum, we kept pairs with relative momentum less 
then 60 MeV/c when computing the source function [24j| . 
We investigated whether this choice of affects the source 
function, but found no noticeable dependence. Because 
of the low temperature of the single particle source, we 
include all simulated pairs (925260 altogether) in the cal- 
culation. In Fig. 12 we see that the resulting source func- 
tion is both non-Gaussian and non-spherical as shown in 
the outward and longitudinal curves. 

Our test correlation is shown in Fig. [31 in the black 
points. We have chosen a much finer binning than is 
usually done in correlation analyses to gain resolution in 
the expanded correlation as we describe in the following 
section. One can clearly see that this test data has not 
been Coulomb corrected in any way and that the cor- 
relation is not spherically symmetric. Furthermore, one 
can see that the tail in the outward direction has pushed 
most of the correlation in the outwards direction into the 
Coulomb hole. 



THE IMAGING PROCEDURE 

In an experiment with resolution Aq = 10 MeV/c, 
there are 40 bins in each direction if one measures to rel- 
ative momentum difference q ~ 200 MeV/c. This means 
that there are 40 3 =64 x 10 3 elements in the data set. If 
we use a Cartesian grid for representing the source with 
10 bins in each direction, then we will need to perform 
64 x 10 3 x 10 3 = 64 x 10 6 six-dimensional integrals just 
to create the inversion matrix. Clearly it is to our benefit 
to work to reduce the size of the problem. Fortunately, 
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FIG. 3: The input model correlation compared with the cor- 
relation expanded in spherical harmonics, stopping at order 
imax = 2. In the panels, we show 5 MeV wide slices cut along 
each axis. 



the full three-dimensional imaging problem can be bro- 
ken up into a series of Id inversions which are much easier 
to handle numerically. 



In this section, we first describe how we convert 
the full three-dimensional problem into a series of one- 
dimensional inversion problems by expanding the corre- 
lation and source in spherical harmonics and the kernel in 
Legendre polynomials. Next, we explain how we expand 
the data itself in this basis. Following this, we explain our 
procedure for setting cut-off value (7 max for all the radial 
integrals and the cut-off order £ max in our expansions. 
Fourth, we detail the representation of the source func- 
tion and we go into some detail how the kernel determines 
the location of the knots in the basis spline expansion of 
the source. The knot locations follow from the Sampling 
Theorem in Fourier theory and require some discussion. 
Finally, we describe the imaging itself. 



Converting the 3d problem into a series of Id 
problems 

We begin by expanding the source and correlation in 
spherical harmonics: 



(10) 



=o m=-e 



and 



ft(q) = V4^g YtmWKtmiq)- (11) 

1=0 m=-£ 

We also expand the kernel in Legendre polynomials in 
/i = cos(6»q !r ): 



A(q,r) = J2P e (ti)K e (q,r). 



(12) 



After inserting these into the Koonin-Pratt equation, we 
simplify the result with the Spherical Harmonic Addition 
Theorem. Doing so, the Koonin-Pratt equation becomes 
a system of linearly independent integral equations: 



Kl m {q)=4* drr z K e (q,r)S em {r). (13) 
Jo 

Instead of attempting a full three-dimensional inversion, 
we will attack this system of Id inversions using the imag- 
ing techniques developed in Refs. 0,0,0- 

Since we focus on like- meson pairs in this paper, we 
write the kernels that we need. For noninteracting spin- 
bosons the kernel K\(q,r) is 



(-l) A/2 jA(2gr) for A even 







otherwise. 



(14) 



If the bosons do interact, the kernel becomes a bit more 
complicated: 

K x (q, r)=J2 ye(r)y*Ar) ( q o ) ~ 6x0 (15) 

We typically truncate the sum over Ps at several times 
feaxfmax/Sc, where q max and r max are the maximum val- 
ues we will need for a particular data set. Here, the object 
in the parentheses is a Wigncr 3 — j symbol and the yi 
are the symmetrized solutions of the radial Klein-Gordon 
equation: 



( V2(-l)^ 2 (2l+l)e i,7t Fg(ri;p)/p even £ 
~ \ odd £ 



(16) 

where at is the Coulomb phase, J) (77; p) is the Coulomb 
radial wave as defined by Messiah |26j , p = qr /he and 77 = 
a QED ^/<7 2 + m 2 /q. While we must include the Coulomb 
interaction between the pair, we may neglect the short- 
range meson exchange potential. 



Converting the 3d data into Id data 

Experimentally measured correlations are inevitably 
binned in q and usually in fixed size Cartesian bins. In 
other words, rather than measuring a continuous corre- 
lation, lZ(q), one measures the average value of the cor- 
relation over a bin: 



n 



ii) 



i 



v. 



d\K{ci). 



(17) 



Here the bin index (i) is shorthand for the set of three 
indices {*Sj*0 5 *l} and the volume of the (i) th bin as 
V(i) . For fixed sized binning, Vu\ is constant and given 
by V(i) = AqsAqoAqL. If we introduce a "window func- 
tion," w^(q), where 



we can write 



«>(i)(q) 



1 if q in Vi 



otherwise, 



(18) 



(19) 



For like pairs, we can reduce the number of data points 
by a factor of two if we notice that the correlation data 
must be invariant under particle label switching. In rel- 
ative coordinates, this is equivalent to: 



(20) 



If we align the edges of the bins along the coordinate 
axes, then we can reflect half of the data back on itself, 
simultaneously reducing the size of the data set by half 
and doubling the statistics on the other half. For the 
purpose of this paper, we will assume all q-bins with 
qs < have been properly reflected into the qs > bins. 

We can now convert a binned correlation in Cartesian 
coordinates to a series of binned correlations in spherical 
coordinates. We invert Eq. l|llf) and average over radial 
bins: 



1 rqj+Aqj/2 [• 

ft<mi=A-/ dq dfiqy;j?)K(q). (21) 

£yc lj Joi — An,- 12 J4tt 



rqj+Aqj/2 
'ij J qj -A qj /2 

Inserting this into l|19|) . wc find 



K 



£mj 



(22) 



(0 



where the transform matrix M is 

I rn+Aij/ 2 r 
Mi mj (i) = — dq drjqF/ m (g)w w (q). 

^1] Jqj-Aqj/2 J Air 

(23) 

Now, the data is a collection of binned one-dimensional 
correlations suitable for imaging. This procedure has the 
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FIG. 4: Terms in the spherical harmonic expansion of the 
model correlation. Since this is a like pair correlation, all odd 
m terms are zero. Because of the symmetry built into our 
model correlation, all terms have only real parts. 



side effect of introducing radial bin-to-bin correlations 
since the Cartesian bins are "non-local" in spherical co- 
ordinates. In fact, the uncertainties are non-trivially dis- 
tributed among all the Im's since the bins in the data are 
re-used to compute each term. Maintaining these corre- 
lations would require us to treat all terms simultaneously, 
resulting in a numerically intractable problem. To make 
progress, we ignore the cross-^m correlations, but main- 
tain the bin-to-bin correlations. They are quantified in 
the covariance matrix of the radial correlation: 

/A 2 H£mj]> = '^^Mi m j(i ) Mt m f(i ) (/A'R(i ) ) 2 . (24) 

W 

In other words, the uncertainties in the correlation data 
are non-trivially distributed among the radial bins in the 
expanded correlation. 

We apply this scheme to the test correlation. In Fig. 0] 
we show the individual terms in the expansion. We only 
show i < 4 terms, although we generated all terms up 
to I = 10 from the test data. In Fig. 0] we recognize 
two important features. First, we see that the Coulomb 
hole is only noticeable in the im = 00 term. The higher 
I terms also have a dip at the origin but this is a result 
of the analyticity of each term of the correlation rather 
than a Coulomb induced effect. The second thing we no- 



tice is that the statistical fluctuations in the expanded 
data arc greatly reduced. This results from the extensive 
angular averaging that occurs in the low £ terms. Thus, 
by picking a seemingly unreasonable binning in q in the 
original correlation, we can achieve a much higher resolu- 
tion in the individual terms of the expanded correlation. 
As we go up in t, the fluctuations become more intense 
as we begin to resolve the angular structure of the start- 
ing bins. This will require us to cut off our spherical 
harmonic expansion at a relatively low £. 

In Fig. we show the (reassembled) expanded correla- 
tion in the squares alongside the original in the circles. It 
is clear that truncating the spherical harmonic expansion 
at a relatively low I (2 in this case) has not distorted the 
correlation and has resulted in significant smoothing of 
the correlation. Some of the smoothing may actually be 
removing information from the correlation, but for the 
most part the smoothing removes angular noise. 

We must truncate our source expansion at some angu- 
lar scale £ m ax- Care must be taken in choosing this £ max : 
if W is too low, then the higher angular scale structure 
in the source will be lost. Actually, the problem is worse 
than this - this higher angular structure can be mistaken 
for noise in the data l28t , interfering with the an- 
gular components we otherwise seek. On the other hand, 
if W is too large, other problems may arise. A large 
4m may be large enough to sample the high frequency 
angular structure imposed by the Cartesian binning at 
large q. Also, since the number of pairs often drop at q 
increases, statistical fluctuations in the high-q bins may 
result in angular aliasing. For a discussion of aliasing, 
see either Ref. ^| or [2£|. In practice, we can set both 
£ max and g max simply by visually inspecting the terms 
in the correlation expansion. In the next subsection we 
will introduce a coarser scheme for setting £ max and g max 
that is nevertheless interesting because of its relation to 
the space-averaged phase-space density. 

We conclude this subsection by describing the physical 
meaning of the terms in the expansion: 

£ — 0: Angle integrated shape. The £m = 00 term is the 
qw correlation if one is analyzing a like pair corre- 
lation. This term is sensitive mainly to Ri nv . 

£ = 1: Dipole distortion. The longitudinal and outwards 
directions are sensitive to any emission time offset 
(i.e. the "Lednicky offset" 31]). If the temporal 
offset occurs in the longitudinally co-moving frame, 
the offset will show up in the £m =11 term. If 
the offset occurs in the colliding nuclei's frame, this 
offset will show up in both the £m =10 and lm = 
11 terms. These terms are absent for like pairs. 

£ = 2: Quadrupole distortion. These terms give the main 
shape of the correlation so are sensitive to the de- 
viation of i? s ide, ^out, and i?iong from i?i nv - The 
lm = 20 term gives i?i on g while linear combina- 
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FIG. 5: "Bump volume" of the model correlation. 



tions of the £m = 00, 20, and 22 terms give i? s idc 
and i? out . 

£ — 3: Octapolc or triaxial distortion. These terms are 
responsible for making the source look something 
like a "boomerang" and are absent for like pairs. 

£ = 4: Hexadecapole distortion. These terms are respon- 
sible for making the source more "box-like." 



Integrals of the data 

In this section, we introduce two integral measures of 
the correlation data: the "bump volume" and the spec- 
tral power. These measures can be used to estimate 
both £ max and g max , but are perhaps more useful as a 
gross characterization of the structure of the full three- 
dimensional correlation. 

The "bump volume" is simply the volume integral of 
the correlation: 



V(q max ) = 47r 



dqq 2 K 00 {q) 



(25) 



This corresponds to the "bump volume" |32j and is di- 
rectly related to the space averaged phase-space density 
[l3j . (/(p)), for identical non-interacting pairs through 



(/(P)> = ± 



Ed 3 N. 



(2s + l)m d 3 p 



-V(<h 



(26) 



This relation is a direct consequence of Eq. (17) in Ref. 

F3| and the sum rule of the correlation described in Ref. 

33j. Unfortunately, as noted in Ref. Q, this sum ride 
can not be applied reliably to data in the presence of 
final-state interactions. Nevertheless, we have plotted 
this volume in Fig. [5] for the test problem. We note that 
the integral is essentially flat past q = 35 MeV/c, indicat- 
ing there is no information accessible in the correlation 
past this momentum. In practice we may have to individ- 
ually tune g max for each term in the spherical harmonic 
expansion of the correlation. 
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FIG. 6: Angular power spectrum of the model correlation 
function. 



A related measure of the correlation is the spectral 
power: 



Vt(q max ) = 4-7T 



E 



dqq 2 \1Zi m {q)Y 



(27) 



This integral gives an indication as to how much of the 
correlation is in each I in the Yi m expansion of the corre- 
lation. The correlation is squared in Eq. (|27f) . allowing 
us to assess the angular content of a correlation even if it 
dips below 1. In Fig. El we show the power spectrum of 
the model correlation. It flattens past t = 4 indicating 
there is no angular information past this order. In fact, 
for this test correlation, it is safe to truncate at an even 
more aggressive £ = 2. 



varying over some range, it makes sense to have very few 
knots there. 

We would like to engineer a basis spline expansion 
such that the approximated source function agrees with 
the true source function at a set of collocation points 
x = {r , ri, r„_i} on the interval [r min , r max ]. In other 
words, 



S{r k ) = ^5 ?; B l (r fc ) 



(29) 



In general, finding the optimal knots t = 
{to, ti, tjVfe-i} that minimize the distance between 
the original function and the approximation in some 
function space metric is a surprisingly difficult problem. 
De Boor |35j presents several schemes for determining 
these optimal knots and all of these schemes begin with 
the nearly optimal choice of 



'i-k-i 



for < i < k - 1 
— for k < i < n 
for n + 1 < i < n 



1 



(30) 

Because De Boor states that this choice is often close 
enough to the optimal knots for practical purposes, we 
adopt these knots for our work. 

Our problem is somewhat more complicated in that we 
want the optimal knots for the convolution of the fitted 
function with some other function, namely the kernel of 
the Koonin-Pratt equation. In the next subsection, wc 
will describe how we use the Sampling Theorem to find 
the collocation points. 



Representing the Source 

As in Ref. 01 > we write the radial dependence in of 
the source functions in terms of basis splines (denoted by 
Bj) 0: 

S tm {r) = Y^S jtm Bj{r) (28) 
i=i 

Basis splines are piece-wise polynomials which are com- 
monly used to approximate functions. Basis splines arc a 
type of spline that has some computational benefits over 
more commonly used splines such as cubic splines |35j . 
A detailed description of basis splines is contained in 

mm. 

To a large extent, the quality of any spline approxima- 
tion of a function is determined by the locations of the 
knots. Knots are the points where the polynomials that 
make up each spline are pasted together. If there is some 
sharp feature in the function to be approximated, then 
we need a high knot density there to resolve the feature. 
On the other hand, if the function is smooth and slowly 



How the kernel determines the collocation points 

The choice of the collocation points for the source and 
bins for the correlation play a surprisingly important role 
in setting the quality of the imaging results. This issue 
was explored in Ref. [T3 . flol ] leading to the development 
of the "optimized discretization" technique. This tech- 
nique can be costly to apply, so we investigate a com- 
putationally cheaper alternative approach. What may 
be surprising is that the choice of binning in the data 
also plays a large role in shaping the quality of the in- 
version. This problem is known in astronomical imaging 
HI IM El Hi and has its root in the Sampling Theo- 
rem in Fourier theory. In this section, we will outline the 
conclusions for the binning from the Sampling Theorem, 
but we will defer a detailed study to a future publication. 

In practice, both TZi rn and Si m are bandwidth lim- 
ited meaning that they are effectively zero outside some 
range. For 7Ze m , this is easy to see. After some high q, 
namely cfecaie, the data is fluctuating about zero and is 
statistically indistinguishable from zero. Similarly, Se m 
is cut-off at some high r sca i e because most particles are 
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emitted from the reaction zone with finite extent. For 
those particles that are produced from resonance decays, 
the decay rate drops exponentially and at some point the 
production rate is effectively zero. 

These bandwidth limits have important implications 
that we now work through. Since each term in the cor- 
relation expansion has compact support, we may expand 
them in terms of spherical Bessel functions, with the re- 
quirement that 7?.£ m (g sca ie) = 0: 



m3t 



(/scale 



(31) 



n=0 

Here, ae n is the n th zero of the I th spherical Bessel func- 
tion. Using the Koonin-Pratt equation we can determine 
the coefficients in this expansion: 
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Here r) = g/g SC aic- 

We can see the meaning of this equation by considering 
the simple case of non- interacting spin-0 bosons: 



1 



*?scalc (j£+l(aen)) 2 



Si m {rin), (33) 



where re n — ai n /2q sca \ c . In other words, the value of the 
source at a few specific points completely determines the 
correlation. Since the source also has compact support: 
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Using this, we can show that 
{-If 12 



S, 



Im 



air, 



2r : 



scale 



(34) 



(35) 



In other words, the source itself is completely determined 
by the correlation at a few points. Both Eqs. I|33|l and 
(|35|l are statements of the Sampling Theorem. 
Why are the points 



and 



air. 
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so important? Well, if we specialize to the I = m = 
case, we see that the size bins in the g-space are related to 
''scale (and the the size bins in the r-space are related to 
(/scale) in a simple manner owing to the fact that oyo = jir: 



Aq = 7r/2r sca i 
Ar = 7r/2<7 sca io 



(38) 



In other words, the largest scale in one space determines 
the smallest resolved scale in the Fourier transform space. 
Eqs. l|3l)|l - l|3*%|) pose an interesting chicken and egg prob- 
lem; these equations imply that there is an optimal bin- 
ning for a correlation given its underlying source. 

Note that the collocation points will be equally spaced 
for £ = 0, but not for higher £. In fact, at low r, the higher 
£ kernels simply do not have the resolution to determine 
the source: all spherical bessel functions go to zero as 
r^Olike j e (r) =r l j£\\ j^. This means that there is a 
"dead zone" at r — for all £ except £ = 0. This dead 
zone occurs in the general case as well. This tells us that 
the collocation points should be more spaced out at low 
r. 

In practice, interactions distort the kernels from the 
simple Fourier kernels. Since the zeros of the kernel de- 
termine the spacing of the collocation points in Eqs. I|36() 
and H37fl for simple Fourier kernels, we expect that the ze- 
ros of more complicated like-pair kernels to play the same 
role. For the like-meson kernel with Coulomb turned on, 
the zeros obtain a q dependence and this is easily taken 
into account. For non-identical pairs, there are no zeros 
in the wavefunctions and our method will fail. We will 
investigate alternative strategies for non-identical pairs 
at a later date. 



3D Imaging 

By combining Eqs. (|17J) and (|13J) . we convert the 
Koonin-Pratt equation into a matrix equation: 

71 = K- S. (39) 

We proceed as in Refs. 0,0,0,113 an d apply Bayesian 
imaging. In Bayesian imaging, we use Bayes Theorem to 
represent the probability of a particular source represent- 
ing the correlation data j3^|- This probability distribu- 

2 

tion is Gaussian, i.e. oc e~ x so the source that maximizes 
the probability minimizes the x 2 '- 

x 2 = (K ■ s - nf ■ (A 2 ^)- 1 • (K ■ S — 11) 



(40) 



The source that docs this is: 



(41) 



S = A 2 S-K T ■ (A 2 ^- 1 K 

The covariance matrix of this source is: 

A 2 S*=(A" T -(A 2 ^)" 1 -A')- 1 . (42) 

In order to stabilize the inversion, we can take advan- 
tage equality constraints |39l | . An equality constraint is a 
condition on the vector of source coefficients that has the 
generic form C • S = c. One example of such a constraint 
is that the source has zero slope at the origin, in which 
case the £m = 00 component obeys 



N M 



SUr^O)=J2 SooiB i( r ^°) =0 - 



(43) 
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With this, the elements of C are Coch — B[(r — > 0) and 
c has one component given by coo = 0. A larger list of 
possible equality constraints and their physical origin is 
tabulated in Tabled 

Equality constraints correspond to a specific form of 
prior information that are easily included in Bayes The- 
orem [l5L l38l| . Adding this prior information amounts to 
adding a penalty term to the % 2; 



X 



A(C ■ S - c) 5 



(44) 



Here A is a trade-off parameter and we may vary it in 
order to emphasize stability in the inversion (by making 
A huge) or to emphasize goodness-of-fit (by setting A to 
zero) . Such an ability to trade-off stability for goodness- 
of-fit is discussed in Numerical Recipes |3(| in detail. 
With this modification of the \ 2 , the imaged source is 

S = A 2 S-(if T -(A 2 ft)- 1 -ft + AC T -c), (45) 

and the covariance matrix of source is 

A 2 S = (K T • (A 2 ^)- 1 • K + XC T -Cy 1 . (46) 

An alternative approach is to use Lagrange multipliers to 
force the constraints to be obeyed. We have investigated 
this approach and found the results to be equivalent. 

IMAGING ANALYSIS OF TEST PROBLEM 
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FIG. 7: Angle averaged imaged and original sources, com- 
pared to best-fit Gaussian source. On the left we plot them 
on a linear scale and on the right on a logarithmic scale. 
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We now turn to the analysis of the test problem. We 
will do this in three stages. In the first, we will analyze 
the £m = 00 term in detail since we can perform cross- 
checks of this term with the CorAL code |4fJ and because 
this term contains several interesting physics elements. 
In the second stage, we will summarize the analysis of 
the higher I terms. In the final stage, we collect all of 
the terms and present the final results of this test. 

The im = 00 term 

The £m = 00 term is the angle averaged, q[ nv , correla- 
tion. Here, we describe the imaging results and compare 
it to Gaussian fits performed with CorAL. We compare 
the restored and the fitted sources to the input model 
source in Fig and to the input model correlation in 
Fig |S1 In the model source, one can clearly see a non- 
Gaussian tail caused by the time dependence of the halo 
term. 

We produced the imaging results shown in Figs. 
and [S] using knots set using the Sampling Theorem with 
<Zmax = 39 MeV/c and 8 3 rd degree basis spline co- 
efficients. We cut the image off at 40 fm and used 
the constraints in Table [I] to control the source behav- 
ior at the origin and at 40 fm. We obtained a final 



FIG. 8: Angle averaged restored correlation, original correla- 
tion and correlation from best-fit Gaussian source. 

X 2 /d.o.f. = 186/39. In Fig. we see that the Gaussian 
core is well resolved and we can begin to resolve the non- 
Gaussian tail at least an order of magnitude down from 
the core. If we convert the source height and half width to 
equivalent Gaussian numbers, we find A = 0.704 ± 0.160 
and i?i IW = 5.77 ± 0.37 fm. Finally, we note the excel- 
lent agreement between the original correlation and the 
correlation produced by uninverting the imaged source. 

We also used CorAL code ^3 to directly fit a Gaussian 
source to the model correlation. CorAL works by folding 
a Gaussian source with the full pion kernel and varying 
the Gaussian source parameters to maximize agreement 
with correlation data. This approach is markedly better 
than the more commonly used Coulomb corrections used 
in correlation analysis. The Gaussian source that CorAL 
uses is: 

s < r > = <v^^-ii;)' (47) 

This parameterization is chosen to conform to radii in 
emission function in Eq. JSJ and is consistent with stan- 
dard radius definitions (i.e. for Id Gaussian tt° correla- 
tion, C(Q inv ) = 1 + Acxp(-Q 2 nv i? 2 „v))- We find best fit 
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Constraint 


Purpose 


Functional form 


r — is a maximum of 5*(r) 
(for like pairs only) 


Constrain the higher £ components 
that are not well controlled due to 
the r dependence of terms in the 
spherical harmonic expansion. 


dS " m (r -0) -0 W,m 
Or 

S tm {r -> 0) = W,m, I 


S(r) = at r = r max 


Smooth oscillations in the source at 
high r caused by aliasing of 
statistical and experimental noise 
in correlation. 


SWn(r max ) =0 W, m 


S(r) — is flat as r — > r max 


Smooth oscillations in the source at 
high r caused by aliasing of 
statistical and experimental noise 
in correlation. 


dg <m (fmax) =0 W, m 



TABLE I: Equality constraints on the Basis Spline representation for non-spherically symmetric sources. Constraints on higher 
order derivatives are possible, but we use 3 rd degree splines so cannot apply these higher order derivative constraints. 



values of i? inv = 5.380 ± 0.134 and A = 0.647 ± 0.029, 
consistent with image. We comment that the Gaussian 
fit to the correlation is slightly worse than that of the im- 
age as we see in both the high q fall off and the Coulomb 
hole in Fig [HI 

Comparing both source extraction methods, we see 
that both the fit and the image reproduce the radii of the 
input source, but both are approximately 20% low at the 
lowest r. In both cases, the sensitivity of of r 2 Ke(q, r) in 
Eq. (|13|) reduces our ability to extract information about 
the source as r — > 0. Improving the input correlation 
with higher statistics improves the agreement between 
input and extracted sources but is not always realistic. 
Making the data binning smaller can improve agreement 
with data but the statistics of the data may not support 
a finer binning. 

Given the excellent agreement between the Gaussian 
fits and the source images, we wondered it if is possible 
to extract the fraction of the source in the core directly 
from the fits or from the equivalent Gaussian parameters 
from the image. If the Gaussian is representative of the 
core in the single particle source, we would expect that 
Aw/ 2 in Eq. ©. Given the fitted A, we find / = 0.804± 
0.018. The fact that / comes out 25% high is due to 
non-Gaussianness of the core itself. The i? ou t and i?i on g 
directions are both strongly distorted because of the time 
emission in the model and the Gaussian fits are not as 
good as they appear. The source integral is a better 
measure of the source except that it cannot separate out 
tail and core. Integrating, we find that the integral is 
0.658 ± 0.009 out to 25 fin and 0.743 ± 0.015 out to 40 
fm, meaning 65.8% of pairs are emitted less than 25 fm 
apart and 74.3% or pairs are emitted less than 40 fm 
apart. 



The higher I terms 

We have imaged the source coefficients up to I — 4. In 
all cases, we used the Sampling Theorem knots, tuned the 
choices of r max and (fecaic to minimize the x 2 to the data, 
and used all the constraints in Table [I] Our parameters 
are summarized in Table ITT1 

Together these terms do not modify the r = fm be- 
havior of the full 3d source because of the leading r 2+i 
behavior of kernel and integration measure. They do 
modify the large r fall-offs at larger r and hence the 
source radii. We show the restored correlations in Fig.El 
They are not exceptionally enlightening but do show the 
overall good agreement below q ma x- We do not use data 
points with q > <7 max , so any agreement there is fortu- 
itous but not required. However, we note the rather high 
X 2 and poor fit to the £m = 22 term. We attribute this 
to the non-Gaussian in the tail in the outwards direc- 
tion causing the relatively sharp rise in this term of the 
correlation, which we have difficulty resolving. We have 
not expanded the test source in lf m 's so do not show the 
source term comparison. 

Collecting the terms: full 3d correlation and source 

We now collect the source and correlation terms and 
compare them to the full 3d input data. We show the cor- 
relation in Fig^J The restored correlation is consistent 
with input, but is significantly smoother. The agreement 
between them is excellent for the higher £ terms and not 
bad for the 00 term. We reproduce both the Coulomb 
hole and the large q fall-off. 

Finally we turn the imaged source shown in Fig 1111 
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FIG. 9: Comparison of the input and restored correlation terms with £ > 0. The input correlation is shown in black symbols 
and the restored source in the open symbols. 
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TABLE II: Parameters used in the imaging analysis of the 
test problem. 



Looking first in the sidewards direction, we clearly see 
the Gaussian core alone. In the outwards and longitudi- 
nal directions, we also clearly see the non-Gaussian tails 
caused by the time dependence of the halo term in Eq JHJ). 
As in the <7i nv correlation, we have the same 20% short 
fall at low r because of the £, m = 00 term. Reading off 
the half widths and source height, we can compute equiv- 
alent Gaussian radii. We find R out = 7.02 ±0.25, i? s ide = 
3.42 ± 0.21, i? long = 5.04 ± 0.31 and A = 0.442 ± 0.067. 
Interestingly, A surprisingly close to the f 2 we put in to 
the test source. The radii are a bit higher than the radii 
we put into the single particle source in Eq. ©. This is 
a result of the mixing of the space and time directions in 



the construction of the source function in Eq. Q ■ In the 
end, it is clear that we not only can image reliably in 3d, 
but we can do it with the same precision as we could in 
one dimension. 

CONCLUSION 

In this paper, we have introduced an extension to the 
source imaging method from Refs. [l3l IlH ll5T| for imag- 
ing full three-dimensional correlations. This extension 
relics on spherical harmonic expansion of the correla- 
tion and source functions, effectively converting an in- 
tractable three-dimensional inversion to a set of simple 
one-dimensional inversions. We have demonstrated this 
technique on a seemingly simple model in which a non- 
Gaussian and non-spherical source tail is generated from 
a simple Gaussian single particle source with finite life- 
time. Within this model, we have shown that we can 
resolve the core radii to a level comparable to Gaussian 
fits. We have also shown that we go far beyond this by re- 
solving the detailed aspherical shape of any non-Gaussian 
tails of the source. 

While we have shown the techniques efficacy only for 
like-pion correlations, imaging is a general tool and is 
limited only by the availability of data and our ability 
to compute relative wavefunctions. Imaging provides a 
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FIG. 10: The input model correlation compared with the cor- 
relation expanded in spherical harmonics, stopping at order 

^max — 4. 
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FIG. 11: The model source function constructed from the 
model emission functions in Eq. The imaged source is 

the curve with the grey error band. We produced the imaged 
source using only the r — constraints and we truncated the 
expansion at i max — 4. 



means to measure source shapes and sizes for complicated 
final states where a Gaussian fit is impossible, such as for 
non-identical correlations. Finally, this extension to the 
imaging technique should shed light on the non-Gaussian 
tails now being seen in advanced analyses of RHIC two- 
particle correlations [4l| . 
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